Aggregrate Multiply-genotyped indivs and assign rockfish with rubias

To start off with, let’s load data and libs:

library(tidyverse)

Loading tidyverse: ggplot2
Loading tidyverse: tibble
Loading tidyverse: tidyr
Loading tidyverse: readr
Loading tidyverse: purrr
Loading tidyverse: dplyr
Conflicts with tidy packages ----------------------------------------------------------------
filter(): dplyr, stats
lag():    dplyr, stats

library(CKMRsim)
library(stringr)
library(readxl)
#meta <- readRDS("../data/processed/meta-data-tibble.rds") # not worrying about meta data right now...
genos <- readRDS("../data/processed/called_genos_na_explicit.rds") #%>%
  #filter(NMFS_DNA_ID %in% meta$NMFS_DNA_ID)  # drop those we don't have meta data for
samples <- readRDS("../data/processed/sample-sheet-tibble.rds") #%>%
  #filter(NMFS_DNA_ID %in% meta$NMFS_DNA_ID)

Some initial filters

Take highest read-depth call for multiply-genotyped DNA_IDs

I’m not sure if there are any of these, but best to leave it in here… particularly for the re-genotyped fish from Brittany’s OSU samples.

Now, here is a harder operation: if an individual is multiply-genotyped, take the genotype with the highest total read depth.

# slow-ish function to get the total read depth column
tdepth <- function(a, d) {
  if(any(is.na(a))) {
    return(NA)
  }
  if(a[1]==a[2]) {
    return(d[1])
  } else {
    return(d[1] + d[2])
  }
  
}
# this takes the highest read-depth instance of each duplicately-genotyped individual.
geno_one_each <- genos %>%
  group_by(NMFS_DNA_ID, locus, gtseq_run, id) %>%
  mutate(total_depth = tdepth(allele, depth)) %>%
  ungroup() %>%
  arrange(NMFS_DNA_ID, locus, total_depth, gtseq_run, id, depth) %>%
  group_by(NMFS_DNA_ID, locus) %>%
  mutate(rank = 1:n()) %>%
  ungroup() %>%
  filter(rank <= 2)

Remove the 6 loci which Hayley has been removing

# read in a list of the 6 loci
to_remove <- read_csv("../data/loci_to_remove.csv")

Parsed with column specification:
cols(
  locus = col_character()
)

# only keep the loci that are not those 6
keepers <- geno_one_each %>%
  anti_join(., to_remove, by = "locus")
# that should leave 90 loci  

Toss out indivs with missing data at more than 15 loci

Now, toss out any individual with fewer than 75 non-missing loci

no_hi_missers <- keepers %>% 
  group_by(NMFS_DNA_ID) %>%
  filter(sum(!is.na(allele)) >= (75*2))

So, we started with 1674 and after filtering out indivs with fewer than 75 genotyped loci, we were left with 1655 individuals. Those are the ones that we will run through rubias to identify to species.

Read in baseline genotypes and remove loci and individuals with too much missing data

# get the data frames into the same format
no_hi_missers2 <- no_hi_missers %>%
  dplyr::select(NMFS_DNA_ID, locus, gene_copy, allele) %>%
  rename(indiv = NMFS_DNA_ID) %>%
  mutate(sample_type = "mixture") %>%
  mutate(repunit = NA) %>%
  mutate(collection = "osu_samples")
# reorder
no_hi_missers2[, c(5:7,1:4)]

# combine the data into a single df
alleles <- bind_rows(spp.id2, no_hi_missers2)
alleles

We are going to do this by turning alleles into integers and spreading it and then getting it into the right format to run rubias.

# first make integers of the alleles
alle_idxs <- alleles %>% 
  #dplyr::select(NMFS_DNA_ID, locus, gene_copy, allele) %>%
  group_by(locus) %>%
  mutate(alleidx = as.integer(factor(allele, levels = unique(allele)))) %>%
  ungroup() %>%
  arrange(indiv, locus, alleidx) # rubias can handle NA's, so no need to change them to 0's
  
# select just the columns to retain and spread the alleles
alle_idx2 <- alle_idxs[,-7]
  
two_col <- alle_idx2 %>%
  unite(loc, locus, gene_copy, sep = ".") %>%
  spread(loc, alleidx)
two_col

# write this file to a thing that can be read-into other softwares
#two_col %>%
# write_csv("csv_outputs/genos_two_col.csv")

Okay, after some reading, it looks like I need to use infer_mixture in rubias. Which require two separate data frames, one with the reference genotypes and the other with the mixture.

I’ll split the data frame that I created, but it needed to be bunged together for the conversion of alleles to integers.

# split up the reference and mixture data frames
sp_mix <- two_col %>%
  filter(sample_type == "mixture")
sp_ref <- two_col %>%
  filter(sample_type == "reference")

# Now that the data are in the corret format, load Rubias
library(rubias)
# perform self-assignment of reference samples
ref_self <- self_assign(sp_ref, gen_start_col = 5)

Summary Statistics:

1403 Individuals in Sample

90 Loci: Plate_1_A01_Sat_GW603857_consensus.1, Plate_1_A11_Sat_GE820299_consensus.1, Plate_2_A09_Sat_EW986980_consensus.1, Plate_2_C08_Sat_EW987116_consensus.1, Plate_2_G06_Sat_EW987118_consensus.1, Plate_3_C03_Sat_GE798118_consensus.1, Plate_4_E10_Sat_EW976030_consensus.1, Plate_4_G06_Sat_EW976181_consensus.1, tag_id_1049.1, tag_id_108.1, tag_id_1184.1, tag_id_1229.1, tag_id_1272.1, tag_id_1366.1, tag_id_1428.1, tag_id_143.1, tag_id_1441.1, tag_id_1449.1, tag_id_1471.1, tag_id_1498.1, tag_id_1558.1, tag_id_1576.1, tag_id_1598.1, tag_id_1604.1, tag_id_1613.1, tag_id_162.1, tag_id_1652.1, tag_id_170.1, tag_id_1708.1, tag_id_1748.1, tag_id_1751.1, tag_id_1762.1, tag_id_179.1, tag_id_1804.1, tag_id_1808.1, tag_id_1810.1, tag_id_1836.1, tag_id_1850.1, tag_id_1880.1, tag_id_1889.1, tag_id_1915.1, tag_id_1950.1, tag_id_1961.1, tag_id_1966.1, tag_id_1982.1, tag_id_1994.1, tag_id_1999.1, tag_id_2008.1, tag_id_2009.1, tag_id_2017.1, tag_id_2062.1, tag_id_2082.1, tag_id_2114.1, tag_id_2134.1, tag_id_2155.1, tag_id_2178.1, tag_id_2182.1, tag_id_220.1, tag_id_2203.1, tag_id_221.1, tag_id_2214.1, tag_id_2237.1, tag_id_2247.1, tag_id_2258.1, tag_id_2301.1, tag_id_2319.1, tag_id_2368.1, tag_id_2499.1, tag_id_250.1, tag_id_2607.1, tag_id_2635.1, tag_id_265.1, tag_id_325.1, tag_id_402.1, tag_id_410.1, tag_id_436.1, tag_id_55.1, tag_id_572.1, tag_id_67.1, tag_id_770.1, tag_id_788.1, tag_id_843.1, tag_id_855.1, tag_id_874.1, tag_id_875.1, tag_id_879.1, tag_id_913.1, tag_id_942.1, tag_id_981.1, tag_id_987.1

48 Reporting Units: Salutus, Satrovirens, Sauriculatus, Saurora, Sbabcocki, Scarnatus, Scaurinus, Schlorostictus, Schrysomelas, Sconstellatus, Scrameri, Sdallii, Sdiaconus, Sdiploproa, Selongatus, Semphaeus, Sentomelas, Sflavidus, Sgoodei, Shelvomaculatus, Shopkinsi, Sjordani, Slevis, Smaliger, Smelanops, Smelanostomus, Sminiatus, Smystinus, Snebulosus, Snigrocinctus, Soculatus, Sovalis, Spaucispinis, Spinniger, Sproriger, Srastrelliger, Sreedi, Srosaceus, Sruberrimus, Srubrivinctus, Srufus, Ssaxicola, Ssemicinctus, Sserranoides, Sserriceps, Sumbrosus, Swilsoni, Szacentrus

48 Collections: Salutus, Satrovirens, Sauriculatus, Saurora, Sbabcocki, Scarnatus, Scaurinus, Schlorostictus, Schrysomelas, Sconstellatus, Scrameri, Sdallii, Sdiaconus, Sdiploproa, Selongatus, Semphaeus, Sentomelas, Sflavidus, Sgoodei, Shelvomaculatus, Shopkinsi, Sjordani, Slevis, Smaliger, Smelanops, Smelanostomus, Sminiatus, Smystinus, Snebulosus, Snigrocinctus, Soculatus, Sovalis, Spaucispinis, Spinniger, Sproriger, Srastrelliger, Sreedi, Srosaceus, Sruberrimus, Srubrivinctus, Srufus, Ssaxicola, Ssemicinctus, Sserranoides, Sserriceps, Sumbrosus, Swilsoni, Szacentrus

5.03365803437079% of allelic data identified as missing

# and take a quick look at those assignments
ref_self %>%
  filter(inferred_repunit == repunit) %>%
  filter(scaled_likelihood > 0.95)

# perform mixture-assignment on baseline colony samples
mix_assign <- infer_mixture(reference = sp_ref, mixture = sp_mix, gen_start_col = 5, method = "MCMC", reps = 2000, burn_in = 100)

Collating data; compiling reference allele frequencies, etc.   time: 1.50 seconds
Computing reference locus specific means and variances for computing mixture z-scores   time: 0.30 seconds
Working on mixture collection: osu_samples with 1436 individuals
  calculating log-likelihoods of the mixture individuals.   time: 0.30 seconds
  performing 100 burn-in and 2000 more sweeps of method "MCMC"   time: 1.79 seconds
  tidying output into a tibble.   time: 0.17 seconds

# That was fast, let's take a look
head(mix_assign)

$mixing_proportions

$indiv_posteriors

$mix_prop_traces

$bootstrapped_proportions

# the individual data is in 
mix_assign$indiv_posteriors %>%
  arrange(desc(log_likelihood))

# are there any rockfish that don't seem like the correct species is in the reference?
osu_assignments <- mix_assign$indiv_posteriors %>%
  filter(PofZ == 1)
osu_assign2 <- osu_assignments[,-10]
#write_csv(osu_assign2, "csv_outputs/osu_assignments.csv")
osu_assign2

The posterior means of group membership in each collection is in the PofZ column - there are 1,220 individuals with a PofZ = 1.

Maybe some of these samples just had a PofZ larger/smaller than 1? (eventually I will update brit_data with all of the appropriate meta data for this project.)

# I think there are individual assignments for each sample to each reference.
# I want to take only the top assignment for each sample.
kept_assignments <- mix_assign$indiv_posteriors %>%
  group_by(indiv) %>%
  filter(PofZ > 0.95) %>%
  arrange(desc(PofZ))
kept_assignments

Using Hayley’s threshold of 0.95 for the PofZ, we have 1,435 individuals retained.

Now take a look at which fish from Brittany’s samples are included:

# using the data I read in previously
tmp2 <- kept_assignments %>%
  left_join(., brit_data, by = c("indiv" = "NMFS_DNA_ID"))
brit_assign <- tmp2[1:21]
# which of brittany's samples still don't have genotype info?
brit_data %>%
  anti_join(., kept_assignments, by = c("NMFS_DNA_ID" = "indiv"))

That didn’t increase the number of samples retained by much…

Are there really 375 samples that don’t have genotypes?

all_assign <- mix_assign$indiv_posteriors %>%
  group_by(indiv)
# which of brittany's samples don't have genotype info?
brit_data %>%
  anti_join(., all_assign, by = c("NMFS_DNA_ID" = "indiv"))

Okay, those 375 samples just aren’t in the data set. They’re not just filtered out by the 0.95 PofZ assignment threshold.

head(brit_data) # here is the "SPECIES" column

kept_assignments %>%
  left_join(., brit_data, by = c("indiv" = "NMFS_DNA_ID")) %>%
  select(1:9,11:30) %>%
  write_csv("csv_outputs/osu_kept_assignments.csv")
  

Repository meta data for Brittany’s samples

batch_4792 <- read_csv("../extdata/batch4792.csv")

Missing column names filled in: 'X16' [16], 'X18' [18]Parsed with column specification:
cols(
  .default = col_character(),
  SAMPLE_ID = col_integer(),
  BATCH_ID = col_integer(),
  LENGTH = col_double(),
  WEIGHT = col_double()
)
See spec(...) for full column specifications.
number of columns of result is not a multiple of vector length (arg 1)700 parsing failures.
row # A tibble: 5 x 5 col     row   col   expected     actual                       file expected   <int> <chr>      <chr>      <chr>                      <chr> actual 1    16  <NA> 22 columns 20 columns '../extdata/batch4792.csv' file 2    17  <NA> 22 columns 20 columns '../extdata/batch4792.csv' row 3    18  <NA> 22 columns 20 columns '../extdata/batch4792.csv' col 4    19  <NA> 22 columns 20 columns '../extdata/batch4792.csv' expected 5    20  <NA> 22 columns 20 columns '../extdata/batch4792.csv'
... ................. ... .............................................................. ........ .............................................................. ...... .............................................................. .... .............................................................. ... .............................................................. ... .............................................................. ........ ..............................................................
See problems(...) for more details.

batch_4969 <- read_csv("../extdata/batch4969.csv")

Missing column names filled in: 'X16' [16], 'X18' [18]Parsed with column specification:
cols(
  .default = col_character(),
  BATCH_ID = col_integer(),
  LENGTH = col_double(),
  WEIGHT = col_integer()
)
See spec(...) for full column specifications.
number of columns of result is not a multiple of vector length (arg 1)242 parsing failures.
row # A tibble: 5 x 5 col     row   col   expected     actual                       file expected   <int> <chr>      <chr>      <chr>                      <chr> actual 1    16  <NA> 22 columns 21 columns '../extdata/batch4969.csv' file 2    17  <NA> 22 columns 21 columns '../extdata/batch4969.csv' row 3    18  <NA> 22 columns 21 columns '../extdata/batch4969.csv' col 4    19  <NA> 22 columns 21 columns '../extdata/batch4969.csv' expected 5    20  <NA> 22 columns 21 columns '../extdata/batch4969.csv'
... ................. ... .............................................................. ........ .............................................................. ...... .............................................................. .... .............................................................. ... .............................................................. ... .............................................................. ........ ..............................................................
See problems(...) for more details.

# somehow the structure for the sample IDs is different in the two files...
# change that.
batch_4792$SAMPLE_ID <- as.character(batch_4792$SAMPLE_ID)
# bind those things together
meta <- bind_rows(batch_4969, batch_4792)

The meta file should contain all the samples for Brittany’s project.

192 samples. That’s still an unfortunate lot. I wonder if they are the same samples that didn’t genotype well last time?

Still 192 samples!

Of primary concern is whether we actually genotyped these at all. Let’s go back to the data.

Yup. The 192 samples are in the hi-missers data - they had lots of missing genotypes.

Looking at the success of re-genotyped samples

# just the samples from gtseq 63
gtseq63 <- samples %>%
  filter(gtseq_run == "gtseq63")

# now add assignments just for those samples
gtseq63 %>%
  left_join(., kept_assignments, by = c("NMFS_DNA_ID" = "indiv")) %>%
  select(1:2, 5:7, 17:22) %>%
  #write_csv("csv_outputs/gtseq63_assignments.csv")
  filter(ssBOX_ID == "R374") #%>%
  #write_csv("csv_outputs/gtseq63_R374_assignments.csv")

Here’s what Cassie said was in that run: This run has box R373 and two mixed plates (1 of those mixed plates is partial). Plate 1 has the samples from plate R374 that need species ID – Brittany’s samples in 1A - 3A; samples for Helen Killeen in 3B - 4A; two re-picks from your NFS project in 4B and 4C; and ONE tentative Mexican rockfish in 4D. It also has re-runs for Brittany’s project. The two NSF re-picks are highlighted in the attached list of black and yellow sibs that you gave me. We re-extracted these since there was tissue remaining and we also re-ran the original DNA extracts in this run. Plate 2rr are all re-runs including your NSF black and yellow sibs and the rest of the failed samples from Bailey’s run for Brittany’s project.

gtseq63 %>%
  filter(ssBOX_ID == "R282")

# which samples from gtseq63 don't have kept assignments?
gtseq63 %>%
  anti_join(., kept_assignments, by = c("NMFS_DNA_ID" = "indiv")) %>%
  filter(ssBOX_ID == "R133")

All the samples from R373 were not re-runs. (23 samples) The 4 samples that failed from R374 were for Helen Killeen.

47 -27 = 20 samples that failed from the re-runs? Was that out of 108 total?

And if we did re-run all of Brittany’s fish from previous runs, is the majority of the missing data from the new samples?